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Abstract. The present paper introduces an efficient and accurate numerical scheme for the solution of a highly 
anisotropic elliptic equation, the anisotropy direction being given by a variable vector field. This scheme is based on 
an asymptotic preserving reformulation of the original system, permitting an accurate resolution independently of the 
anisotropy strength and without the need of a mesh adapted to this anisotropy. The counterpart of this original procedure 
is the larger system size, enlarged by adding auxiliary variables and Lagrange multipliers. This Asymptotic-Preserving 
method generalizes the method investigated in a previous paper 1111 to the case of an arbitrary anisotropy direction field. 

1. Introduction 

Anisotropic problems are common in mathematical modeling of physical problems. They occur in 
various fields of applications such as flows in porous media [3l I20j , semiconductor modeling |29j , quasi- 
neutral plasma simulations [10], image processing [37, 38 , atmospheric or oceanic flows }36| and so on, 
the list being not exhaustive. The initial motivation for this work is closely related to magnetized plasma 
simulations such as atmospheric plasma |24[ 126] , internal fusion plasma [51 112) or plasma thrusters [T] . 
In this context, the media is structured by the magnetic field, which may be strong in some regions and 
weak in others. Indeed, the gyration of the charged particles around magnetic field lines dominates the 
motion in the plane perpendicular to magnetic field. This explains the large number of collisions in the 
perpendicular plane while the motion along the field lines is rather undisturbed. As a consequence the 
mobility of particles in different directions differs by many orders of magnitude. This ratio can be as 
huge as 10 10 . On the other hand, when the magnetic field is weak the anisotropy is much smaller. As 
the regions with weak and strong magnetic field can coexist in the same computational domain, one 
needs a numerical scheme which gives accurate results for a large range of anisotropy strengths. The 
relevant boundary conditions in many fields of application are periodic (for instance in simulations of 
the tokamak plasmas on a torus) or Neumann boundary conditions (atmospheric plasma for example 
[S])- For these reasons we propose a strongly anisotropic model problem for wich we wish to introduce 
an efficient and accurate numerical scheme. This model problem reads 

{-V-AV0 e = / in Q, 
n-AV0 e = O ondfijv, (1.1) 
cj) e = Q ondil D , 

where fl C M 2 or f2 C K 3 is a bounded domain with boundary <90 = dflo U dfljy and outward normal n. 
The direction of the anisotropy is defined by a vector field B, where we suppose divB = and B^Q. 
The direction of B is given by a vector field b = B/\B\. The anisotropy matrix is then defined as 

A = -A\\b®b+(Id-b®b)Ax(Id-b®b) (1.2) 

£ 

and <9f2 j = {i£ <9i7 | b(x) -n = 0}. The scalar field A\\ >0 and the symmetric positive definite matrix 
field A±_ are of order one while the parameter < e < 1 can be very small, provoking thus the high 
anisotropy of the problem. This work extends the results of [TTJ , where the special case of a vector field 
b, aligned with the z-axis, was studied. An extension of this approach is proposed in [6] to handle more 
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realistic anisotropy topologies. It relies on the introduction of a curvilinear coordinate system with one 
coordinate aligned with the anisotropy direction. Adapted coordinates are widely used in the framework 
of plasma simulation (see for instance [H [21 [31]), coordinate systems being either developped to fit 
particular magnetic field geometry or plasma equilibrium (Euler potentials 35 , toroidal and poloidal 
[Ulllll, quasiballooning [15], Hamada [19] and Boozer [7] coordinates). Note that the study of certain 
plasma regions in a tokamak have motivated the use of non-orthogonal coordinates systems |21) . In 
contrast with all these methods, we propose here a numerical scheme that uses coordinates and meshes 
independent of the anisotropy direction, like in [53] • This feature offers the capability to easily treat time 
evolving anisotropy directions. This is very important in the context of tokamak plasma simulation, 
the anisotropy being driven by the magnetic field which is time dependent. 

One of the difficulties associated with the numerical solution of problem (|1.1|) lies in the fact that 
this problem becomes very ill-conditioned for small 0<e^;l. Indeed, replacing e by zero yields an 
ill-posed problem as it has an infinite number of solutions (any function constant along the b field 
solves the problem with e = 0). In the discrete case the problem translates into a linear system which 
is ill-conditioned, as it mixes the terms of different orders of magnitude for £<gl. As a consequence 
the numerical algorithm for solving this linear system gives unacceptable errors (in the case of direct 
solvers) or fails to converge in a reasonable time (in the case of iterative methods). 

This difficulty arises when the boundary conditions supplied to the dominant 0(l/e) operator 
lead to an ill-posed problem. This is the case for Neumann boundary conditions imposed on the part 
of the boundary with b-n^O as well as for periodic boundary conditions. If instead, the boundary 
conditions are such that the dominant operator gives a well-posed problem, the numerical difficulty 
vanishes. One can resort to standard methods, as the dominant operator is sufficient to determine the 
limit solution. This is the case for Dirichlet and Robin boundary conditions. The problem addressed 
in this paper arises therefore only with specific boundary conditions. It has however a considerable 
impact in numerous physical problems concerning plasmas, geophysical flows, plates and shells as an 
example. In this paper, we will focus on Neumann boundary condition since they represent a larger 
range of physical applications. The periodic boundary conditions can be addressed in a very similar 
way. 

Numerical methods for anisotropic problems have been extensively studied in the literature. Dis- 
tinct methods have been developed. Domain decomposition techniques using multiple coarse grid 
corrections are adapted to the anisotropic equations in [XT] [27] . Multigrid methods have been studied 
in |16[ 132] . For anisotropy aligned with one or two directions, point or plane smoothers are shown to 
be very efficient [28] . The hp- finite element method is also known to give good results for singular per- 
turbation problems [30) . All these methods have in common that they try to discretize the anisotropic 
PDE as it is written and then to apply purely numerical tricks to circumvent the problems related to 
lack of accuracy of the discrete solution or to the slow convergence of iterative algorithms. This leads 
to methods which are rather difficult to implement. 

The approach that we pursue in this paper is entirely different: we reformulate first the original 
PDE in such a way that the resulting problem can be efficiently and accurately discretized by straight- 
forward and easily implementable numerical methods for any anisotropy strength. Our scheme is related 
to the Asymptotic Preserving method introduced in |22) . These techniques are designed to give a precise 
solution in the various regimes with no restrictions on the computational meshes and with additional 
property of converging to the limit solution when e — > 0. The derivation of the Asymptotic Preserving 
method requires identification of the limit model. In the case of Singular Perturbation problems, the 
original problem is reformulated in such a way that the obtained set of equations contain both the limit 
model and the original problem with a continuous transition between them, according to the values 
of e. This reformulated system of equation sets the foundation of the AP-scheme. These Asymptotic 
Preserving techniques have been explored in previous studies, for instance quasi-neutral or gyro-fluid 
limits [HI [H], as well as anisotropic elliptic problems of the form (jl.ip with vector b aligned with a 
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coordinate axis [TTJ [B] . 

In this paper, we present a new algorithm which extends the results of The originality of 

this algorithm consists in the fact, that it is applicable for variable anisotropy directions b, without 
additional work. The discretization mesh has not to be adapted to the field direction 6, but is simply a 
Cartesian grid, whose mesh-size is governed by the desired accuracy, independently on the anisotropy 
strength e. All this is possible by a well-adapted mathematical framework (optimally chosen spaces, 
introduction of Lagrange multipliers). The key idea, as in is to decompose the solution <f) into 
two parts: a mean part p which is constant along the field lines and the fluctuation part q consisting 
of a correction to the mean part needed to recover the full solution. Both parts p and q are solutions 
to well-posed problems for any e > 0. In the limit of e — > the AP-reformulation reduces to the so 
called Limit model (L- model) , whose solution is an acceptable approximation of the P-model solution 
for E-Cl (see Theorem I2.2[) . In [IT] the Asymptotic Preserving reformulation of the original problem 
was obtained in two steps. Firstly, the original problem was integrated along the field lines (z-axis) 
leading to an e-independent elliptic problem for the mean part p. Secondly, the mean equation was 
subtracted from the original problem and the e-dependent elliptic problem for the fluctuating part 
q was obtained. This approach however is not applicable if the field b is arbitrary. In this paper we 
present a new approach. Instead of integrating the original problem along the arbitrary field lines, we 
choose to force the mean part p to lie in the Hilbert space of functions constant along the field lines 
and the fluctuating part q to be orthogonal (in L 2 sense) to this space. This is done by a Lagrange 
multiplier technique and requires introduction of additional variables thus enlarging the linear system 
to be solved. This method allows to treat the arbitrary b field case, regardless of the field topology 
and thus eliminates the limitations of the algorithm presented in We note that an alternative 

method, bypassing the need in Lagrange multipliers, is proposed in [8|. It is based on a reformulation 
of the original problem as a fourth order equation. 

The outline of this paper is the following. Section [2] introduces the original anisotropic elliptic 
problem. The original problem will be referred to as the Singular- Perturbation model (P-model). The 
mathematical framework is introduced and the Asymptotic Preserving reformulation (AP-model) is 
then derived. Section |3] is devoted to the numerical implementation of the AP-formulation. Numerical 
results are presented for 2D and 3D test cases, for constant and variable fields b. Three methods are 
compared (AP-formulation, P-model and L-model) according to their precision for different values of e. 
The rigorous numerical analysis of this new algorithm will be the subject of a forthcoming publication. 

2. Problem definition 

We consider a two or three dimensional anisotropic problem, given on a sufficiently smooth, bounded 
domain f2cR d , g?=2,3 with boundary d£l. The direction of the anisotropy is defined by the vector 
field be(C°°(fl)) d , satisfying |6(a:)| = l for aUatefi. 

Given this vector field b, one can decompose now vectors ueR d , gradients V</>, with (j>(x) a scalar 
function, and divergences V-i>, with v(x) a vector field, into a part parallel to the anisotropy direction 
and a part perpendicular to it. These parts arc defined as follows: 

vu := (v-b)b, v±:=(Id — b(gib)v, such that v = v\\ +vj_ , 

V||fl!>:= {b-V<j>)b, V±0:= (Id-b®b)S7(f), such that V0 = V||0 + V_l0, (2.1) 
Vii -v := V'Uii , Vj_ • v := V -v± , such that V-?; = V| -v + V± -v, 

where we denoted by <g> the vector tensor product. With these notations we can now introduce the 
mathematical problem, the so-called Singular Perturbation problem, whose numerical solution is the 
main concern of this paper. 
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2.1. The Singular Perturbation problem (P-model) We consider the following Singular 
Perturbation problem 

r -IV r (A||V||<^)-Vx-(AxVx^) = / in ft, 

(P) I ±n v (A ll \7 ll (j) e )+n ± -(A ± V ± (l) e ) = on dfl N , (2.2) 
{(f) e = ondn D , 
where n is the outward normal to ft and the boundaries are defined by 

dn D = {xedn | b(x)-n = 0}, dn N = dQ\dn D . (2.3) 

The parameter < e < 1 can be very small and is responsible for the high anisotropy of the problem. 
The aim is to introduce a numerical scheme, whose computational costs (simulation time and memory) , 
for fixed precision, arc independent of e. 

We shall assume in the rest of this paper the following hypothesis on the diffusion coefficients and the 
source terms 

o 

Hypothesis A Let /G£ 2 (ft) and dVtr>^0. The diffusion coefficients A^EL°°(fl) and A± G 
^dxd{L°°(Q)) are supposed to satisfy 

0<A <A\\(x)<Ai, fa.a.xett, (2.4) 
A)IM| 2 <^_l(2><^iIM| 2 , \/v£R d andf. a. a. x G ft. (2-5) 



As we intend to use the finite element method for the numerical solution of the P-problem, let us put 
(|2.2[) under variational form. For this let V be the Hilbert space 

V:={^Gif 1 (^)/^ D =0}, (0,^) v :=(V||0,V|^)^+£(Vx</),VxV)^- 

Thus, we are seaking for <jf G V, the solution of 

a\\(^^)+ea ± ( ( /> e ^) = s(f,iP), W>GV, (2.6) 

where (•,•) stands for the standard L 2 inner product and the continuous bilinear forms an :Vx V->R 
and a± : V x V — > R are given by 

a||(0,V>) := / A||V||^-V||^<fa;, ax(&V0 : = / (A±V ± <P)-V ± ^dx. (2.7) 
Jn Jn 

Thanks to Hypothesis A and the Lax-Milgram theorem, problem (12.21) admits a unique solution <fi £ G V 
for all fixed e > 0. 

2.2. The Limit problem (L-model) The direct numerical solution of (|2.2I) may be very 
inaccurate for eel. Indeed, when e tends to zero, the system reduces to 

{-V||-(A||V||0)=O in ft, 

ji||.(A||V||<A)=0 on<9ftjv, (2.8) 

(f) = on 9ft_o. 
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This is an ill-posed problem as it has an infinite number of solutions 4"EG, where 

G = {<t>eV | V||^ = 0}, (2.9) 

is the Hilbcrt space of functions, which are constant along the field lines of b. This shows that the 
condition number of the system obtained by discretizing (|2.2[) tends to oo as e — > so that its solution 
will suffer from round-off errors. 

For this reason, we should approximate (|2.2[) in the limit e — > differently. Supposing that (ff — > f> Q 
as e — > we identify (at first formally) the problem satisfied by <fP . From the above arguments we know 
that (f>° £G- Taking now test functions ip € Q in (|2.6j) . we obtain 

/ A ± V ± <f) e ■V ± ipdx= [ fipdx. (2.10) 
Jn Jo. 

Passing to the limit e —> into this equation yields the variational formulation of the problem satisfied 
by cj) (Limit problem): find (j>° EG, the solution of 

(L) [ A ± V ± <t) -\7 ± ipdx= [ fipdx ,VipeG, (2.11) 
Jo. Jo. 

which is a well posed problem. Indeed, the space G C V is a Hilbert space, associated with the inner 
product 

(0,^)0 :=(Vx^VxV)l», V^gS, (2.12) 
and the norm || • \ \g is equivalent to the H 1 norm. This is due to the Poincare inequality, as 
ll^lli 2 <C||V0||i 2 = C||V||0|| 2 L2 + C||V ± 0|| 2 L2 = C||Vx0||! 2 , WeG- 

Hypothesis A and the Lax-Milgram lemma imply the existence and uniqueness of a solution <fi° € G of 
the Limit problem (|2.1ip . 



Remark 2.1. Let us restrict ourselves for the moment to the simple special case (considered in a pre- 
vious paper Ulf ) of the two dimensional domain Q = (0,L x ) x (0,L Z ) in the (x,z) plane with a constant 
b-field aligned with the Z-axis: 

b=(Tj- (2.13) 

The functions in the space G are independent of z so that G can be identified to Hq(0,L x ). The limit 
problem (|2.1ip now reads: Find (jP in H^{0,L x ) verifying 

A ± (x)d x( j) {x)d x i;{x)dx= / f{x)^{x)dx, ^eH^(0,L x ), 



where Aj_(x) = (1/ L z ) J^ z A± t n(x,z)dz and f(x) = (1/L Z ) f(x,z)dz are the mean values of Aj_ and 
f along the field lines. The limit solution </>° thus verifies a one- dimensional elliptic equation whose 
coefficients are integrated along the anisotropy direction: 



-d x (A ± (x)d x <f>°(x)^=f(x) on (0,1,), 
0°(O) = 0°(L ;c ) = O. 
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(2.14) 



We see now that <j> (x) is a solution to the one dimensional elliptic problem so that it belongs to 
H 2 (0,L x ) provided f^L 2 (Q). Since (fr as a function of (x,z) does not depend on z, we have also 
(f>° €H 2 (Sl). This conclusion (<fp £ H 2 (Q)) remains valid in the case of a cylindrical three dimensional 
domain £1 = Q xy x (0,L Z ) in the (x,y,z) space with any sufficiently smooth Q xy in the (x,y) plane and the 
field b aligned with the Z -axis, b= (0,0, 1)* . Indeed, it is easy to see that <p =<f) (x,y) solves in this case 
an elliptic two dimensional problem in tt xy similar to \2.1J$ so that we can apply the standard regularity 
results for the elliptic problems. These examples show that it is reasonable to suppose <fP GH 2 (fl) also 
in more general geometries of Q and b. This can be indeed proved under the hypotheses in Appendix 
HI by specifying the (d— 1) dimensional elliptic problem for </>° . The proof being rather lengthy and 
technical, we prefer to postpone it to a forthcoming work \13]j . 

2.3. The Asymptotic Preserving approach (AP-model) In this section we introduce the 
AP-formulation, which is a reformulation of the Singular Perturbation problem (|2.2I) . permitting a 
"continuous" transition from the (P)-problem (|2 .2[) to the (L)-problem (I2.1ip , as e — > 0. For this purpose, 
each function is decomposed into its mean part along the anisotropy direction (lying in the subspace Q 
of V) and a fluctuating part (cf. [IT]) lying in the L 2 -orthogonal complement A of Q in V, defined by 

A:={4>GV |(M0 = O, V^ee?}. (2.15) 

Note that (•,•) denote here and elsewhere the inner product of L 2 (il). 
In what follows, we need the following 

Hypothesis B The Hilbert- space V admits the decomposition 

V = G® ± A ) (2.16) 

with Q given by V2. 9\) and A given by V2.15\) and where the orthogonality of the direct sum is taken with 
respect to the L 2 -norm. Denoting by P the orthogonal projection on Q with respect to the L 2 inner 
product: 

P:V^Q such that {P(j) 1 ^) = {(j),^) V0eV,V>e£, (2.17) 
we shall suppose that this mapping is continuous and that we have the Poincare-Wirtinger inequality 

||0--Pfc ( n)<C||V||0|| L 2 ( n), V^eV. (2.18) 

Applying the projection P to a function <f> is nothing but a weighted average of <f> along the anisotropy 
field lines of b. The space Q is the space of averaged functions (the parallel (gradient of these averaged 
functions being equal to zero), whereas the space A is the space of the fluctuations (the Average of the 
fluctuations being equal to zero). Note that the decomposition (|2.16l) is not self evident and it may in 
fact fail on some "pathological" domains fi. Indeed, although one can always define an L 2 -orthogonal 
projection P<f> on the space of functions constant along each field line, for any <j) with square-integrable 
V||</>, one cannot assure in general that P(f> belongs to V for V since one may lose control of the 
perpendicular part of the gradient of Pip. Fortunately however, Hypothesis B is typically satisfied for 
the domains of practical interest. The interested reader is referred to Appendix [X] for an example of a set 
of assumptions on fi and b which entail Hypothesis B and which resume essentially to the requirement 
for the field b to intersect d£l n in a uniformly non-tangential manner and for the boundary components 
8Hn and dfto to be sufficiently smooth. 
Let us also define the operator 



Q:V^A, Q = I-P. 
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(2.19) 



Each function ip G V can be decomposed uniquely as <fi = p + q, where p = P<p £ Q and q = Q4> G A. Using 
this decomposition, we reformulate the Singular- Perturbation problem (|2.2I) . Indeed, replacing <f> £ :— 
p e + q e in problem (|2.2[) and taking test functions 77 G 5 and ^GA leads to an asymptotic preserving 
formulation of the original problem: Find (p £ ,q £ ) (zG * A such that 

f a ± (p £ ,rj) + a ± (q £ , r?) = (/, 77) Vt7G£, 

^ (2.20) 

i a||(^0+ea^£)+£a ± (p £ ,0 = e(/,0 V£eA 

Contrary to the Singular Perturbation problem (|2.2p , setting formally e = in (12.20[) yields the system 

(a ± { P °,r,) + a x {q\n) = (/,*?), Vt^ 

i (2.21) 

Un(?°,0 =0, VfeA, 

which has a unique solution (p°,q°) & Q x A, where p° is the unique solution of the L-problem p. lip and 
q° = 0. Indeed, taking £ = <Z° as test function in the second equation of (|2.21[) yields V| | g° = 0, which 
means q°GG- But at the same time, q°GA, so that q° <^Q DA — {0}. Setting then q° = in the first 
equation of (|2.21[) . shows that p° is the unique solution of the L-problem. 

Theorem 2.2. For every e>0 the Asymptotic Preserving formulation i2.20\) . under Hypotheses A 
and B, admits a unique solution [p £ ,q £ ) <EG x A, where (f> £ :=p £ + q £ is the unique solution in V of the 
Singular Perturbation model \2.2\) . 
These solutions satisfy the bounds 

W\\m { n)<C\\f\\ LH n), M\\m { n) < C\\f\\ L 2 {n) , \\p £ \\m { a) <C\\f\\ L 2 (a) , (2.22) 

with an e-independent constant C>0. Moreover, we have 

(f> e ^(f> , p £ ^(/> and q £ ^0 in H 1 (ft) as e -> , (2.23) 

where (jp E G is the unique solution of the Limit model h2.11\) . 

Proof. The existence and uniqueness of a solution for the P-problem as well as L-problem are 
consequences of the Lax-Milgram theorem. The existence and uniqueness of a solution of (|2.20l) is then 
immediate by construction, remarking that the decomposition 4> £ —p £ + q £ is unique. 
The bound ||0 e ||//i(n) ^C||/||l 2 (o) is obtained by a standard elliptic argument. Furthermore, p £ = P(jf 
where P is the L 2 -orthogonal projector on G, which is a bounded operator in V by (|1.4j) . This implies 
the estimates for p £ and q £ in (|2.22j) . Since p £ e Q and q £ G„4 are bounded, there exist subsequences 
p £n and q £n that weakly converge for e„ -> to some p° € Q and q° G A. Taking e — £ n in (|2.20p and 
passing to the limit e n — >0 we identify (p°,q°) with the unique solution of (I2.2ip . i.e. p° = <fi° is the 
unique solution of (|2 . 1 1[) and q° = 0. Since the limit does not depend on the choice of the subsequence, 
we have the weak convergence as e — > 0, i.e. 

p £ ^ £ ^oP° in H\n), q £ ^ in If 1 (ft). 

We shall prove now that these convergences are actually strong. Introducing e £ =p £ — p , we have 

a ± {e £ ,n)+a ± {q £ ,r 1 ) = 0, MneG- 

Taking now 7/ = e £ in this relation and adding it to the second equation in (I2.20p . where we put £ = q £ /e, 
yields 

-a n (q £ ,q £ ) + a ± (q £ + e £ ,q £ + e £ ) = (f,q £ )-a ± (p°,q £ ). (2.24) 

e 



Due to the Poincare-Wirtinger equation (I2.18[) . there exist a constant C>0 such that 



Ulm^^Ca^q) 1 / 2 , VqeA. (2.25) 

In combination with a Young inequality this gives (/,<? e ) < ||/||L 2 ll<7 e ||L 2 < £ Tf 1 1 /1 11,2 + | (g e ,g e ). Us- 
ing this in the right hand side of (|2.24l) . we arrive at 

Noting that q £ + e £ — <p £ —p° and V|| e e = we can rewrite this last inequality as 

^a { £ -p ,^ -p a ) + a ± (^ -p ^ -p )<e^\\f\\l 2 -a ± (p ,q £ ). 

Since a±(p°,q £ ) — >Q as e — >0 (thanks to the weak convergence q £ — ^0) we observe that cf> £ — >p° strongly 
in i? 1 (f2). Reminding again that p £ — Pip £ and P is bounded in the norm of iJ 1 (f2), we obtain also 
p £ — > Pp° =p°, which entails q £ — > 0. □ 

Remark 2.3. Let us return to the simple special case discussed in remark [2. 1\ i.e. f2 = (0,L X ) x (0,L Z ) 
and the b-field given by &2.13\) . Remind that the space Q can be identified in this case with the space 
of functions constant along the Z-axis, which means Q '■={(/) G V / d z cj) = 0}. The space A is orthogonal 
(with respect to the L 2 -norm) to Q and thus contains the functions that have zero mean value along 
the Z-axis, i.e. A:—{(j)&V/ J Q * <fi(x,z)dz = 0} . Therefore, for <j) £ =p £ + q £ G V, the function p £ is the 
mean value of <p £ in the direction of the field b: 

p s = ±- [ ~<f> £ dz, (2.26) 

L z Jo 

and q £ is the fluctuating part with zero mean value: 

q £ = (t> £ -^- ( ' <j) £ dz. (2.27) 

Hypothesis B is thus easily verified. The results obtained in this special case were presented in a previous 
paper In the case of an arbitrary b-field, formula i2.26\) is generalized as il.2\) in Appendix [21 

where the length element along the b-field line is weighted by the infinitesimal cross-sectional area of 
the field tube around the considered b-field-line. This formula can be thus interpreted as a consequence 
of the co-area formula. Note that in the special case of a uniform anisotropy direction, the limit 
problem can easily be formulated as an elliptic problem depending on the only transverse coordinates 
(see equation (|2.14p ). The size of the problem is thus significantly smaller than that of the initial one. 
This feature still occurs for non-uniform b-fields as long as adapted coordinates and meshes are used. 
In our case, aligned and transverse coordinates are not at our disposal and the solution of the limit 
problem must be searched as a function of the whole set of coordinates. 

2.4. Lagrange multiplier space The objective of this work is the numerical solution of system 
(|2.20|) and the comparison of the obtained results with those obtained by directly solving the original 
problem (|2.2j) . In a general case, when the field b is not necessarily constant, the discretization of 
the subspaces Q and A, is not straightforward, as in the simpler case [TT]. In order to overcome this 
difficulty a Lagrange multiplier technique will be used. 
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2.4.1. The A space To avoid the use of the constrained space A, we can remark that A can 
be characterized as being the orthogonal complement (in the L 2 sense) of the Q-sp&ce. Thus, instead 
of (I2.20p . the slightly changed system will be solved: find (p £ ,q £ ,l £ ) €G xV x Q such that 



a,, (<f , + ea ± (q £ , £) + ea ± (p £ , £) + (l £ , £) = e(f, V£ e V, 
(g e ,x) = v x e£. 



(2.28) 



The constraint (q £ ,x) — 0, Vx^G 1S forcing the solution q £ to belong to A, and this property is carried 
over to the limit e—tO. We have thus circumvented the difficulty of discrctizing A by introducing a 
new variable and enlarging the linear system. 

PROPOSITION 2.4. Problems &2.20\) and H2.28\) are equivalent. Indeed, (p £ ,q £ ) eQ x A is the unique 
solution of i2.20\) if and only if (jf ,q £ ,l e ) G £7 x V x Q with l £ = is the unique solution of \2.28]) . 

Proof. Let (p £ ,q £ )^G x A be the unique solution of (|2.20[) . Then, it is immediate to show that 
(p £ ,q £ ,0) solves (|2~25|) . Let now (p £ ,q £ ,l £ ) e G x V x Q be a solution of (|2~28|) . Then, the last equation 
of (12.281) implies that q £ € A Choosing in the second equation as test function £ 6 one gets 

eai(a £ ,0 +«u(p £ ,0 + (Z e ,0 = e(/,0 , G 5, 

which because of the first equation in (|2.28p . yields (Z e ,£) = for all £,£G- Thus / £ =0. □ 

2.4.2. The Q space In order to eliminate the problems that arise when dealing with the 
discretization of Q, the Lagrange multiplier method will again be used. First note that 



p£Go 



| V M p = ( J A ll \7 ll p-\7 ll \dx = a ll (p,\) = 0,\/\e£ 

\p^ v IpgV, 



(2.29) 



where £ is a functional space that should be chosen large enough so that one could find for any p G V a 
Ag£ with Vi|A = Viip. On the other hand, the space C should be not too large in order to ensure the 
uniqueness of the Lagrange multipliers in the unconstrained system. A space that satisfies these two 
requirements under some quite general assumptions to be detailed later, can be defined as 



£:={\eL 2 (n) /Vn\eL 2 (n), A, an( „=0}, with dtt in ;={xedn / b(x)-n<0}. 



(2.30) 



Using the characterization (|2.29[) of the constrained space Q, we shall now reformulate the system 
(|2~28l) as follows: Find (p £ , X £ , g E ,( £ ,/i £ )eVx£xVxVx£ such that 

fa ± (p £ ,77)+a ± (q e ,7y) + a||(77,A £ ) = (/,77), Vr^eV, 
a| | (p e ,k)=0, Vkg£, 
(I/ 7 ) { a ll (q £ 7 0+ea ± (q £ ,0+za±(p £ ,0 + (l £ ,0=z(f,0, V£gV, (2.31) 
(q £ ,x)+a\\(x,^) = 0, V X eV, 
I a||(Z e ,r) = 0, WeC. 

The advantage of the above formulation, as compared to (|2.20p . is that we only have to discretize the 
spaces V and C (at the price of the introduction of three additional variables) , which is much easier than 
the discretization of the constrained spaces Q and A. More importantly, the dual formulation (|2.31[) 
does not require any change of coordinates to express the fact that p £ is constant along the 6-field lines 
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and that q e averages to zero along these lines. Therefore this formulation is particularly well adapted 
to time-dependent 6-fields, as it does not require any operation which would have to be reinitiated as 
b evolves. The system (|2.3ip will be called the Asymptotic-Preserving formulation in the sequel. 
To analyse this Asymptotic-Preserving formulation, we need the following 

Hypothesis B' The trace \\dn in is well defined for any AG V as an element of L 2 (dfli n ), with contin- 
uous dependence of the trace norm in L 2 (dfli n ) on ||A||y. Moreover, the Hilbert space 

V = {0£L 2 (r!)/V||0eL 2 (r!)}, (^)9:=(^) + (V||&V||V), (2-32) 
admits the decomposition 

V = g®£, (2.33) 

where Q is given by 

0:={0eV/V||0 = O}, (2.34) 

and C is given by H2. 30\) . The spaces Q and G — QC\V are related in the following way: if g&Q is such 
that J dn r/gda — for all n EQ , then .9 = 0. 

The decomposition (|2.33|) is quite natural. It tells simply that any function (f> can be decomposed on 
each field line as a sum of a function that vanishes at one given point on this line and a constant (which 
is therefore the value of 4> a t this point). Hypothesis B' will be thus normally satisfied in cases of 
practical interest. For example, we prove in Appendix [A] that the set of assumptions on the domain Q, 
and the 6-field which can be used to verify Hypothesis B, is also sufficient (but far from necessary) for 
Hypothesis B'. We are now able to show the relation between systems (|2.28[) and (|2.31j) . 
Proposition 2.5. Assuming Hypotheses A, B and B\ problem \2. 31\) admits a unique solution 
(p e , A e , q E ', l e ', [i e ) GVx£xVxVx£, where (p E ,q e ,l e ) G Q x V x Q is the unique solution of h2.28\) . 

The proof of Proposition 12.51 is based on the following two lemmas 
Lemma 2.6. Assume Hypothesis B' and letp^V be such that a\\(p,X) =0, VAg£. Then p^Q . 

Proof. Take any tjGV and decompose r/ = A+g with XeC and g£G- We have a^(p,g) — 0, hence 
a ||(P; 7 ?) = f° r an ??GV. This entails V||p = 0, hencepG^. □ 

Lemma 2.7. Assume Hypothesis B' and let F GV* be such that F (t)) = for all rj G Q . Then the problem 
of finding A G C such that 

a ll (r,,X)=F(r l ), Vt?gV, (2.35) 

has a unique solution. 

Proof. Consider the bilinear form b on V x V 

b(u,v) = a^(u,v)+ / uvda 

By Hypothesis B', this is an inner product on V. Indeed, if b(u,u) = then uGQnC so that u~0. 
Riesz representation theorem implies that the problem of finding fi^V such that 

b(n,n) = F(r,), Vt?gV, 

has a unique solution. We can now decompose fi = \ + g with A G C and g G £?. This yields 

a \\{ r l^)+ I V9dcr = F(n), VrjeV, 
JdQ in 
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so that, in particular, J dil rjgdo~ = for all 77 G{? which implies <? = 0. We see now that A is a solution 
to (|2.35p . The uniqueness follows easily. □ 
Let us now prove Proposition 12.51 

Proof of existence in Proposition [2751 Take (p £ ,q s ,l £ ) 6?xVx(/ as the unique solution of (|2.28[) . 

Then, equations 2,3,5 in (|2.31[) are immediately satisfied. It remains to choose properly the Lagrange 
multipliers X £ ,/i £ g£ to satisfy equations 1,4 in (12.31[) . For this, let us define Fi,F2 G V* by 

Fi(ri):=±a n (<f,ri), F 2 (r,):=-(q e ,r l ), V77GV. (2.36) 

These functionals are indeed continuous in the norm of V since their definitions do not contain the 
derivatives in directions perpendicular to b. Since F\ (77) = F2 (r?) = for all rjGG, Lemma [2.71 implies 
the existence of X £ G C and /i £ G £, such that 

0||(T ? ,A e ) = F 1 (t J ), ail (x^ £ )=F 2 ( x ), V^xGV. (2.37) 

Taking rj,x G V G V we observe (cf. the second line in (|2.28[) where l e — 0) 

a\\(x,^) = -(q £ ,x), VxGV, 

which coincides with equations 1,4 in (|2.3ip . 

Proof of uniqueness in Proposition ^. 5l Consider the solution to system (|2.31[) with / = 0. Lemma 
[2761 implies then that p e , l E G Q n V = Q and (p s , q e , l E ) G 5 x V x (? verifies ([2728)) with / = so that p e = 
q £ = l e = 0by Proposition [274] Equations 1,4 in (f273T|) now tell us that A £ ,/i £ G Q, but QnC = {0}, hence 
A £ = ^ £ = 0. □ 

The presence of 1 /e in the formulas (|2.36l) , (12.37)) defining X e indicates at a first sight that A 6 may 
tend to 00 as £->0 which would be disastrous for an AP numerical method based on (|2.31[) at very 
small e. Fortunately A e remains bounded uniformly in e in the cases of practical interest. It suffices to 
suppose that the limit solution 0° is in H 2 (Q) which is a reasonable assumption as discussed in Remark 

ED 

Proposition 2.8. Assume Hypotheses A, B, B' and (f)°eH 2 (£l) where (jp is the solution to V2.ll}) . 
Then X e introduced in \2. 31)) satisfies 

||V||A £ || L 2<Cmax(||/|| i2 ,||0°|| H2 ) (2.38) 

with a constant C independent of e. 

Proof. We will denote all the e-independent constants by C in this proof. We start from relation 
(|2.24[) in the proof of Theorem 12.21 Dropping the positive term a±(q e + e £ ,q £ + e £ ) it can be rewritten 
as 

- £ a n (q £ ,q £ )<(f,q £ )-a ± (cb ,q £ ). 

Since </>° Gi? 2 (f2) we can integrate by parts in the integral defining a±(cf) ,q e ): 

- a± (0V) = - / A 1 _V ± <f>°-V ± g e dx 

= - [ {Id-bb t )A 1 V^(jP-nq £ da+ [ (V ± • A ± V^ a )q e dx 
Jdn N Jo. 

<C\\4P\\m{}\q e \\L^dn N) + \\q e \\L H n)) 
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since V<p° has a trace on d£l and its norm in L 2 (9iljv) is bounded by C||0°||#2. Thus, 

i 1 1^11^7- 1 < "^-^1 1 Cq'",^") < C| I ^ 1 1 1 1^7^ 1 1 J^- (f^) H~ l^>° 1 1^- CI 1^" I l^-COf^^) 1 1 l.r.^c^)) - 
By Poincare-Wirtinger inequality (I2.18[) (note that Pq e = 0) and by Hypothesis B' we have 

max(||<f || L 2 (fi) ,||<f || L 2 (a0jv) ) <C||V||(f || L 2 

so that 

i||V||g e |U2<Cmax(||/|| i2! ||/||^). 
This is the same as since V||A £ = ^\\q e according to (|2~M]) and (pTBTj) . □ 

Remark 2.9. The Limit model h2.11\) . reformulated using the Lagrange multiplier technique, now 
reads: Find (0°, A ) E V x C such that 

{[ A 1 _V^cj) -V^dx+ [ A||V||V-V||A°dx= f ftpdx VipEV 
Jn Jn Jn (23g) 

J A ll V ll <f>°-V\\Kdx = Q Vks£. 

Problem 12. 39]) is also well posed assuming Hypotheses A, B, B' and (fp EH 2 (tt). Indeed, the uniqueness 
of the solution to \2. 39}) can be proved in exactly the same manner as in the proof of Proposition \2.5\ 
above. To prove the existence of a solution, it suffices to take the limit e— >0 in the first two lines of 
12.31}) . Indeed, we know by Theorem \ 2.2\ that p £ — >-0° , the solution to i2.ll}) . and q e — >0 in H x (fX). 
Moreover, the family {V||A £ } is bounded in the norm of L 2 (fl) by Proposition ] 2. 8[ We can take therefore 
a weakly convergence subsequence {V||A £ "} and identify its limit with {V||A } with some X°eC (cf. 
Lemma\2~7\) to see that (0° ,A°) eV x C solves 1MB- 

3. Numerical method 

This section concerns the discretization of the Asymptotic Preserving formulation (|2.31[) . based 
on a finite element method, and the detailed study of the obtained numerical results. The numerical 
analysis of the present scheme is investigated in a forthcoming work [13 , in particular we are interested 
in the convergence of the scheme, independently of the parameter e > 0. 

Let us denote by Vh C V and £h<Z £ the finite dimensional approximation spaces, constructed by 
means of appropriate numerical discretizations (see Section [3. II and Appendix B). We are thus looking 
for a discrete solution (p e h , A|, q h , l e h , /i £ h ) € Vh x Ch x Vh xV^x Ch of the following system 

' a±{p h ,r]) + a ± {ql,n) + a ll (i 1 ,X' h ) = (f,r]) , Vr? E V h , 

a\\{p e h ,K) = 0, \fnECh-, 

< a ll (q h ,0 + Ea x (qi,0 + ea x (p h ,0 + (lU)=z(f,0, V£eV ft , (3.1) 

k a||(Z|,r) = 0, VrEC h . 

Our numerical experiments indicate that the spaces Vh and Ch can be always taken of the same 
type and on the same mesh. The only difference between these two finite element spaces lies thus in 
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the incorporation of boundary conditions. In general, let denote the complete finite element space 
(without any restrictions on the boundary) which should be H 1 conforming but otherwise arbitrarily 
chosen. We define then 

V h = {v h eX h /v h \ d n D =0}, (3.2) 



Ch = {Ah <^X h /\ h \gn iri vjdn D =0}. (3.3) 

While this choice of Vh is straight forward, the boundary conditions in Ch require special attention. 
Indeed, nothing in the definition (|2.30l) of space C on the continuous level indicates that its elements 
should vanish on dflo- However, this liberty on drijj is somewhat counter- intuitive. Indeed, the 
Lagrange multiplier X £ S C serves to impose V||p e = for some function p £ taken from the space V. But, 
for p€ V the trace on dfto is zero so that V||p e = there without the help of a Lagrange multiplier. 
Of course, this argument is not valid on the continuous level since the trace of functions in C does not 
even necessarily exist. However, this may become very important on the finite element level. Indeed, 
we provide in Appendix [51 an example of a finite element setting without incorporating Xh\dn D —0 into 
the definition of Ch, which leads to an ill-posed system (|3.1[) . To avoid this difficulty, we choose Ch as 
in (|3.3[) in all our experiments, thus obtaining well-posed problems. 

3.1. Discretization Let us present the discretization in a 2D case, the 3D case being a simple 
generalization. The here considered computational domain f2 is a square f2= [0,1] x [0,1]. All simula- 
tions are performed on structured meshes. Let us introduce the Cartesian, homogeneous grid 

Xi = i/N x , 0<i<iV X) y 3 =j/N y , 0<j<N y , (3.4) 

where N x and N y are positive even constants, corresponding to the number of discretization intervals 
in the x- resp. y-direction. The corresponding mesh-sizes are denoted by h x > resp. h y > 0. Choosing 
a Q2 finite element method (Q2-FEM), based on the following quadratic base functions 



' (x — Xi-2)(x — Xi-\) _ r I 

2hp, X t [Xi—2 ,Xi\, 

(xi-L2—x)(xi + 1 —x) r 1 f) 

^-t X£[Xi,X i+2 \, 1 U Vj 

else 



(y-yj-2){y-yj-i) 



for even i,j and 



x)(x — Xi-i) _ r l 

^± — £2 '- xe[xi-i,x i+ i\ 

else 



(Vj + 2 


2hl 

-y)(yj+i-y) 











-y)(y-yj-i) 








v^[vj,y 3 +2] 

else 



else 



(3.5) 



(3.6) 



for odd we define 



X h ■= {v h =^2vtj 9 Xi (x) 6 Vj (y)} , 



1 -j 



We then search for discrete solutions (p s h , q^, l h ) SV^x Vh x Vh and (A|, fi e h ) € Ch X Ch with Vh and Ch 
defined by (|3.2[) and (|3.3[) . This leads to the inversion of a linear system, the corresponding matrix 
being non-symmetric and given by 



(3.7) 



(A, 


A 


A! 





\ 


A 














sAi 







c 











C 





A 


V 








A 


J 
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£ 


AP scheme 


Limit model 


Singular Perturbation scheme 


L 2 error 


H error 


L 2 error 


H error 


L 2 error 


H error 


1U 


7.2 x 1CT 6 


4.7 x 1CT 3 


5.0 x 10° 


3.51 xlO 1 


7.2 x 10~ 6 


4.7 x 10~ 3 


1 


7.3 x 10" 7 


4.7 x 1CT 4 


5.0 xlO" 1 


3.51x10° 


7.3 x 10" 7 


4.7 x 10" 4 


HT 1 


1.47 x 1(T 7 


9.6x 10~ 5 


5.0 x 10~ 2 


3.51 xlO" 1 


1.45 x 10~ 7 


9.4x 10~ 5 


10~ 4 


1.28 x 10~ 7 


8.3x 10~ 5 


5.0 x 10~ 5 


3.61 xlO~ 4 


1.26 x 10~ 7 


8.2 x 10~ 5 


10~ 6 


1.28 x 1(T 7 


8.3 x 1CT 5 


5.2x 10~ 7 


8.4 x 10" 5 


5.9x 10~ 7 


8.2 x 10~ 5 


10 -io 


1.28 x 10~ 7 


8.3x 10~ 5 


1.28 x 10~ 7 


8.3 x 10" 5 


9.9x 10~ 3 


3.12 xlO" 2 


io- 15 


1.28 x 1CT 7 


8.3x 10~ 5 


1.28 x 10~ 7 


8.3 x 10~ 5 


7.1 x IO" 1 


2.23 x 10° 



Table 3.1 - Comparison between the Asymptotic Preserving scheme, the Limit model and the 
Singular Perturbation model for h = 0.005 (200 mesh points in each direction) and constant 6: 
absolute L 2 -error and i/ 1 -error, for different e-values. 



The sub-matrices A , A\ resp. C correspond to the bilinear forms aii(v), a_i_(v) resp. (•,•), used in 
equations (|2.31[) and belong to R( Ara: + 1 )( JV !'+ 1 ) x ( Ara =+ 1 )( Ar K+ 1 ). The matrix elements are computed using 
the 2D Gauss quadrature formula, with 3 points in the x and y direction: 



J J f{x,y)= UiUjfixi^j), (3.8) 



rl 1 

I 

where = yo = 0, x±\ =y±i = iyf , w o = 8/9 and u)±i = 5/9, which is exact for polynomials of degree 
5. 

3.2. Numerical Results 

3.2.1. 2D test case, uniform and aligned &-field In this section we compare the numerical 
results obtained via the Q2-FEM, by discretizing the Singular Perturbation model (|2.2p , the Limit model 
(|2.1ip and the Asymptotic Preserving reformulation (|2.31[) . In all numerical tests we set Aj_ = Id and 
A\\ = 1. We start with a simple test case, where the analytical solution is known. Let the source term 
/ be given by 

/= (4 + £)7r 2 cos(27rx)sin(7ry)+7r 2 sin(7ry) (3.9) 

and the b field be aligned with the x-axis. Hence, the solution <f> e of (|2.2j) and its decomposition 
4> £ — p s + q £ write 

(j) £ =sm(7ry)+£cos(27nc)sin(7n/), (3.10) 

p e = sm(ny) , q £ = £COs(27ra)sin(7T2/) . (3-H) 



We denote by c/)p, </>l, 4>a the numerical solution of the Singular Perturbation model (|2.2[) . the 
Limit model (|2 . 1 1|) and the Asymptotic Preserving reformulation (|2.31[) respectively. The comparison 
will be done in the L 2 -norm as well as the ii/^-norm. The linear systems obtained after discretization 
of the three methods are solved using the same numerical algorithm — LU decomposition implemented 
in a solver MUMPS [2J. 

In Figure 13.11 we plotted the absolute errors (in the L 2 resp. H ^norms) between the numerical 
solutions obtained with one of the three methods and the exact solution, and this, as a function of 
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(a) L 2 error for a grid with 50 X 50 points. 

10 r ■ 1 1 , 1 ■ 1 1 1 1 ■ 1 1 , 1 ■ 1 - 

<p> — »- 



(c) L 2 error for a grid with 100 X 100 points. 



1e-05 



(c) L 2 error for a grid with 200 X 200 points. 



(b) H 1 error for a grid with 50 X 50 points. 



(P) □ 

ty •••*»• 

(Ap — 9- 




















V 







(d) H error for a grid with 100 X 100 points. 



(P) — it- 
Si 

(AP 



(f) H error for a grid with 200 X 200 points. 



Figure 3.1 - Absolute L 2 (left column) and H 1 (right column) errors between the exact solution 
cf> E and the computed numerical solution 4>a (AP), (j>L (L), 4>p (P) for the test case with constant 
b. The error is plotted as a function of the parameter e and for three different mesh-sizes. 



the parameter e and for several mesh-sizes. In Table 13. 1[ we specified the error values for one fixed 
grid and several e-values. One observes that the Singular Perturbation finite element approximation 
is accurate only for e bigger than some critical value ep, the Limit model gives reliable results for e 
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method 


# rows 


=ff non zero 


time 


X 2 -error 


i? 1 -error 


AP 


50 x 10 3 


1563 x 10 3 


13.212 s 


1.02 x 10~ 6 


3.34 x 10~ 4 


L 


20 x 10 3 


469 x 10 3 


5.227 s 


1.14 xl0~ 6 


3.34 x 10~ 4 


P 


10 x 10 3 


157 xlO 3 


3.707 s 


1.02 xl0~ 6 


3.27 x 10~ 4 



Table 3.2 - Comparison between the Asymptotic Preserving scheme (AP), the Limit model (L) 
and the Singular Perturbation model (P) for h = 0.01 (100 mesh points in each direction) and 
fixed e — 10 -6 : matrix size, number of nonzero elements, average computational time and error 
in L 2 and H 1 norms. 



smaller than el, whereas the AP-scheme is accurate independently on e. The order of convergence for 
all three methods is three in the L 2 -norm and two in the if -norm, which is an optimal result for Q2 
finite elements. When designing a robust numerical method one has therefore two options. The first 
one is to use an Asymptotic Preserving scheme, which is accurate independently on e, but requires the 
solution of a bigger linear system. The second one is to design a coupling strategy that involves the 
solution of the Singular Perturbation formulation and the Limit problem in their respective validity 
domains. This is however a very delicate problem, since we observe that the critical values Ep and el 
are mesh dependent, namely Ep inversely proportional to h and El proportional to h. Therefore for 
small meshes there may exist a range of e-values, where neither the Singular Perturbation nor the Limit 
model finite element approximation give accurate results. For our test case, this is even the case for 
meshes as big as 200 x 200 points, if one regards the L 2 -norm. This mesh-size is generally insufficient 
in the case of real physical applications. 

Another interesting aspect with respect to which the three methods must be compared, is the 
computational time and the size of the matrices involved in the linear systems. Table 13.21 shows that 
the Asymptotic Preserving scheme is expensive in computational time and memory requirements, as 
compared to the other methods. Indeed, the computational time required to solve the problem is almost 
four times bigger than that of the Singular Perturbation scheme. Moreover, the Asymptotic Preserving 
method involves matrices that have five times more rows and ten times more nonzero elements than 
the matrices obtained with the Singular Perturbation approximation. It is however the only scheme 
that provides the /i-convergence regardless of e. In order to reduce the computational costs, a coupling 
strategy for problems with variable e will be proposed in a forthcoming paper. In sub-domains where 
e > Ep the Singular Perturbation problem will be solved, in sub-domains where e < El the Limit problem 
will be solved and only in the remaining part, where neither the Limit nor the Singular Perturbation 
model are valid, the Asymptotic Preserving formulation will be solved. 

3.2.2. 2D test case, non-uniform and non-aligned &-field We now focus our attention on 
the original feature of the here introduced numerical method, namely its ability to treat nonuniform b 
fields. In this section we present numerical simulations performed for a variable field b. 

First, let us construct a numerical test case. Finding an analytical solution for an arbitrary b 
presents a considerable difficulty. We have therefore chosen a different approach. First, we choose a 
limit solution 

4P = s\\\(jty + a{y 2 — y)cos(nx)) , (3.12) 

where a is a numerical constant aimed at controlling the variations of b. For a = 0, the limit solution 
of the previous section is obtained. The limit solution for a = 2 is shown in Figure 13.21 We set a = 2 
in what follows. Since <fi a is a limit solution, it is constant along the b field lines. Therefore we can 
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Figure 3.2 - The limit solution for the test case with variable b. 



determine the b field using the following implication 

V||</>° = => b x - 
which yields for example 



OX 



dy 



a(2y— 1)cos(7tj:) +7T 
Tra(y 2 — y)sin(7rx) 



(3.13) 



(3.14) 



Note that the field B, constructed in this way, satisfies div£? = 0, which is an important property in the 
framework of plasma simulation. Furthermore, we have B ^ in the computational domain. Now, we 
choose cj) £ to be a function that converges, as e— ?>0, to the limit solution (j>°: 



(3.15) 



<ff = sin (ny + a(y 2 —y)cos(nx)) + £CO.s(27ra)sin(7ry) . 
Finally, the force term is calculated, using the equation, i.e. 

/ = -V± • (A x V ± <f ) - ~V|| • V||0 £ ). 

As in the previous section, we shall compare here the numerical solution of the Singular Perturbation 
model (|2.2|) . the Limit model (|2.11j) and the Asymptotic Preserving reformulation (|2.31|) . i.e. <fip, <pL, 
4>A with the exact solution ()3. 15|) . The L 2 and H 1 -errors are reported on Figure 13.31 and Table 13.31 
Once again the Asymptotic Preserving scheme proves to be valid for all values of e, contrary to the 
other schemes. There is however a difference compared to the constant-6 case. For a variable b , the 
threshold value ep seems to be independent on the mesh size and is much larger than that of the 
uniform b test case. This observation limits further the possible choice of coupling strategies, since 
even for coarse meshes there exists a range of e- values, where neither the Singular Perturbation nor the 
Limit model are valid. The coupling strategy, involving all three models, remains however interesting 
to investigate. 
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e 


AP scheme 


Limit model 


Singular Perturbation scheme 


\j error 


H error 


L 2 error 


H"^ error 


L 2 error 


H error 


1 n 

1U 


7.2x 10~ 6 


4.6 x 10~ 3 


5.0 x 10° 


3.50 x 10 1 


7.2 x IO" 6 


4.6 x 10~ 3 


1 


7.1 xlO^ 7 


4.6 x IO" 4 


5.0 x IO" 1 


3.50 x 10° 


7.1 xl0~ 7 


4.6 x IO" 4 


1(T 2 


2.05 x IO" 7 


1.33 xlO~ 4 


5.0 xlO -3 


3.50 x IO" 2 


2.05 xlO" 7 


1.33 XlO" 4 


icr 4 


2.12 x 10" 7 


1.38 xlO~ 4 


5.0x 10~ 5 


3.77 x IO" 4 


1.74 xlO" 6 


1.43 XlO" 4 


10- 7 


2.17 x IO" 7 


1.41 xlO~ 4 


2.22 x IO" 7 


1.41 xlO -4 


1.68 xlO" 3 


1.26 xlO" 2 


io- 10 


2.17 x IO" 7 


1.41 xlO~ 4 


2.17 x IO" 7 


1.41 xlO -4 


3.93 xlO" 1 


1.35 x 10° 


io- 15 


2.17 x IO" 7 


1.41 xlO" 4 


2.17 x IO" 7 


1.41 xlO -4 


6.7 x IO" 1 


2.32x10° 



Table 3.3 - Comparison between the Asymptotic preserving scheme, the Limit model and the 
Singular Perturbation model for h = 0.005 (200 mesh points in each direction) and variable 6: 
absolute L 2 -error and """-error. 



In the next test case we investigate the influence of the variations of the b field on the accuracy of 
the solution. We would like to answer the following question: what is the minimal number of points per 
characteristic length of b variations required to obtain an acceptable solution. For this, let us modify 
the previous test case. Let b = B/\B\, with 

B=( a ±lJ CO ltT2^)^ (3-16) 



mira{y — y)sm(mirx) 
m being an integer. The limit solution and 4> e are chosen to be 

4>° = sm(iry + a(y 2 — y)cos(rrnrx)) , (3-17) 
<j) e = sin (j:y + a(y 2 —y)cos(rmrx)) +ecos(27rx)sin(7ry) . (3.18) 

We perform two tests: first, we fix the mesh size and vary m to find the minimal period of b for which 
the Asymptotic Preserving method yields still acceptable results. We define a result to be acceptable 
when the relative error is less then 0.01. In the second test m remains fixed and the convergence of the 
scheme is studied. The results are presented on Figures 13.41 and 13.51 

For e = l and 400 mesh points in each direction (h = 0.0025) the relative error in the L 2 -norm, 

defined as 4n ^^ l2 '" ) is below 0.01 for all tested values of l<m<50. The relative 7J 1 -error 
-^rr ^iLsllg) exceeds the critical value for m>25. For e = 10- 2 ° the maximal m for which the er- 

N < Mllffl(f2) 

ror is acceptable in both norms is 20. The minimal number of mesh points per period of b variations is 
40 in the worst case, in order to obtain an 1% relative error. 

Figure 13.51 and Table 13.41 show the convergence of the Asymptotic Preserving scheme with respect 
to h for m = 10 and e — IO -10 . We observe that for big values of h the error does not diminish with h. 
Then, for h< 0.025 the scheme converges at a better rate then 2 for if 1 -error and 3 for L 2 -error. For 
h< 0.00625 (160 points) the optimal convergence rate in the iJ 1 -norm is obtained (which is 32 mesh 
points per period of b). The method is super-convergent in the whole tested range for the L 2 -error. 

These results are reassuring, as they prove that the Asymptotic Preserving scheme is precise even 
for strongly varying fields for relatively small mesh sizes, which was not evident. Indeed, the optimal 
convergence rate in the ii^-norni is obtained for 32 mesh points per b period, and an 1% relative error 
for 40 points. It shows that accurate results can be obtained in more complex simulations, such as 
tokamak plasma for example. The application of the method to bigger scale problems is the subject of 
an ongoing work. 
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(a) L 2 error for a grid with 50 X 50 points. 
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(b) H 1 error for a grid with 50 X 50 points. 
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(c) L 2 error for a grid with 100 X 100 points. 



(d) H error for a grid with 100 X 100 points. 



1e-06 - 
1e-07 L 



(c) L 2 error for a grid with 200 X 200 points. 



(f) H error for a grid with 200 X 200 points. 



Figure 3.3 - Absolute L 2 (left column) and H 1 (right column) errors between the exact solution 
(f> E and the computed solution cj>A (AP), 4>l (L), cj>p (P) for the test case with variable b. Plotted 
are the errors as a function of the small parameter e, for three different meshes. 
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1 e-08 I 1 1 1 1 1 1 1 1 1 1 1e-08 I 1 1 1 1 1 1 1 1 1 

5 10 15 20 25 30 35 40 45 50 5 10 15 20 25 30 35 40 45 50 

(a) e = l (b)e = 10- 20 

Figure 3.4 - Relative L 2 and H 1 errors between the exact solution <jf and the computed solution 
4>a (AP) for h = 0.0025 (400 points in each direction) as a function of m and for e — 1 respectively 

io- 20 . 



10 




Figure 3.5 - Relative L 2 and H 1 errors between the exact solution </> e and the computed solution 
4>a (AP) for m = 10 and e = 10~ 10 as a function of h. 




3.2.3. 3D test case, uniform and aligned &-field Finally, we test our method on a simple 
3D test case. Let the field b be aligned with the X-axis: 

(3.19) 

> > 

Let f2= [0,1] x [0,1] x [0,1], and the source term / is such that the solution is given by 

<fi E = sin(7ry)sin(7rz) -|-ecos(27rx)sin(7ry)sin(7rz) , (3.20) 
p e = sin(7r?/)sin(7r2:) , q e = ecos(27ra;)sin(7ry)sin(7rz) . (3-21) 

Numerical simulations were performed on a 30 x 30 x 30 grid. Once again all three methods are com- 
pared. The L 2 and if ^errors are given on Figure |3~51 The numerical results are equivalent with those 
obtained in the 2D test with constant b. Note that it is difficult to perform 3D simulations with more 
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-f+- UU111LO UVjl JJOl 1UU 


L error 


H - error 


U.l 


o 
Z 


a Twin — 1 
4./ X 1U 


1.U5 


U.Uo 


A 

4 


O.Z X 1U 


1 on 


0.025 


8 


1.82 x 10" 1 


4.3 x 10" 1 


0.0125 


16 


1.89 x 10~ 2 


6.4 x 10~ 2 


0.00625 


32 


1.41 xl0~ 3 


1.00 x 10~ 2 


0.0003125 


64 


9.3 x 10~ 5 


2.21 x 10~ 3 


0.0015625 


128 


6.1 x 10~ 6 


5.5 x 10" 4 


0.00078125 


256 


4.6x 10~ 7 


1.36 x 10~ 4 



Table 3.4 - Relative L 2 and H 1 errors between the exact solution (f> £ and the computed solution 
4>a (AP) for m = 10 and e = 10 -10 as a function of h. 



refined grids, due to memory requirements on standard desktop equipment as we are doing now. Every 
row in the matrix constructed for the Singular Perturbation model, can contain up to 125 non zero 
entries (for Q2 finite elements), while matrices associated with the Asymptotic Preserving reformulation 
have rows with up to 375 non zero entries. Furthermore the dimension of the latter is five times bigger. 
The memory requirements of the direct solver used in our simulations grow rapidly. The remedy could 
be to use an iterative solver with suitable preconditioncr. Finding the most efficient method to inverse 
these matrices is however beyond the scope of this paper. In future work we will address this problem 
as well as a parallelization of this method. 




(a) L 2 error for a grid with 30 X 30 X 30 points. (b) H 1 error for a grid with 30 X 30 X 30 points. 



Figure 3.6 - Absolute L 2 (left column) and H 1 (right column) errors between the exact solution 
(jf and the computed solution cj>A (AP), (f>L (L), (f>p (P) for the 3D test case. The errors are 
plotted as a function of the anisotropy ratio e. 



4. Conclusions 

The asymptotic preserving method presented in this paper is shown to be very efficient for the 
solution of highly anisotropic elliptic equations, where the anisotropy direction is given by an arbitrary, 
but smooth vector field b with non-adapted coordinates and meshes. The results presented here gener- 
alize the procedure used in |llj and have the important advantage to permit the use of Cartesian grids, 
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independently on the shape of the anisotropy. Moreover, the scheme is equally accurate, independently 
on the anisotropy strength, avoiding thus the use of coupling methods. The numerical study of this 
AP-scheme shall be investigated in a forthcoming paper, in particular the £-independent convergence 
results shall be stated. 

Another important related work consists in extending our methods to the case of anisotropy ratios 
e, which are variable in f2 from moderate to very small values. This is important, for example, in plasma 
physics simulations as already noted in the introduction. An alternative strategy to the Asymptotic 
Preserving schemes whould be to couple a standard discretization in subregions with moderate e with 
a limit (s— >0) model in subregions with small e as suggested, for example, in [5J 125) . However, the 
limit model is only valid for e <C 1 and cannot be applied for weak anisotropies. Thus, the coupling 
strategy requires existence of a range of anisotropy strength where both methods are valid. This is 
rather undesirable since this range may not exist at all, as illustrated by our results in Fig. 13.11 

Appendix A. Decompositions V — Q®^ A, V = Q®C and related estimates. 

We shall show in this Appendix that all the statements in Hypotheses B and B' can be rigorously 
derived under some assumptions on the domain boundary dfl and on the manner in which it is inter- 
sected by the held b. We assume essentially that b is tangential to dft on dfto and that b penetrates 
the remaining part of the boundary 9f2jv at an angle that stays away from on <9f2jv. We assume also 
that <9fijv consists of two disjoint components for which there exist global and smooth parametrizations. 
This last assumption can be weakened (existence of an atlas of local smooth parametrizations should be 
sufficient) at the expense of lengthening the proofs. The precise set of our assumptions is the following: 

Hypothesis C The boundary of is the union of three components: OVId where b-n = 0, dflin 
where b-n<—a and dVt out where b-n>a with some constant a>0. Moreover, there is a smooth 
system of coordinates on dfli n meaning that there is a bounded domain ri„eM d_1 and 

a one-to-one map hi n 'Ti n —^M. such that hi n <E C 2 (r ,„) and d£li n is the image of /im(£i,---,6i-i) as 
goes over Ti n . The matrix formed by the vectors {dhi n /dS y \,...,dhin/d£ l i,n) is invertible 
for all (£i, . . . G f m . Similar assumptions hold also for dfl ou t (changing f™ to T ou t and hi n to h ou t)- 



Using this hypothesis we can introduce a system of coordinates in such that the held lines of 
b coincide with the coordinate lines. To do this consider the initial value problem for a parametrized 
ordinary differential equation (ODE): 

ay 

aid 

Here is Revalued and £' stands for . For any fixed £' £ r jn , equation should 

be understood as an ODE for a function of £4. Its solution X(£',£d) g° es then over the field line of b 
starting (as £d = 0) at the point on the inflow boundary d£li n , parametrized by This field line hits 
the outflow boundary dU out somewhere. In other words, for any £' eT^ there exists L(£')>0 such 
that X(£',L(£')) €dCl ou t- The domain of definition of X is thus 

D = {{i'X d )eM d I £'er in and <&<£(£')}• 

Gathering the results on parametrized ODEs, from for instance [34], we conclude that X(£',£d) = 
A(£i,...,£d) is a smooth function of all its d parameters, more precisely XgC 2 (D). Evidently, the 
map X is one-to-one from D to fl and thus provide a system of coordinates for J7. Moreover this 

system is not degenerate in the sense that the vectors dX/d^i,...,dX/d^d are linearly independent at 
each point of £1. Indeed, if this was not the case, then there would exist a non trivial linear combination 
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XidX / dt;A \-AddX/d£d that would vanish at some point in f2. But, ODE implies 



f) . d fix d f)x 

so that, the unique solution of this ODE, i.e. the linear combination X4=i^i§fS would vanish on the 
whole field line, in particular on the inflow. But this is impossible since ^ = , i = 1 , . . . , d— 1 on the 
inflow, while f|^ = ^ and the vectors f^g^N ■ ■ ■ , Q^j\ are linearly independent for all € 

Ti„. We see thus that the Jacobian J = det(dX j/d£i) does not vanish on fl so that we can assume that 
m< J < M everywhere on Q with some positive constants m and M (assuming that J is positive does 
not prevent the generality since if J is negative in f2 than one can replace £1 by — £1). Since X e C 2 (f2), 
we have also that JeC 1 (ri). 

One also sees easily that the top of D is given by a smooth function L(£'). Indeed, i s 
determined for each £' G T j n from the equation X(£', L(£') ) = /i out (77) with some 77 = (771 , . . . , r?<j-i ) £ r ou t . 
We know already that this equation is solvable for £ d = L(£'), 771 , . . . , T)d-i for any £' & T in . To conclude 
that the solution depends smoothly on £' we can apply the implicit function theorem to the equation 

F (C ';$d,m> ■ ■ • >%-i) = X (C \€d) ~ h out (m, ■ ■ ■ m-i) = °- 

Indeed, the i? d -valued function F is smooth and the matrix of its partial derivatives with respect to 
£(2,771, •••>%-i is invertible since dF '/ 'd£d = b and dF/di]i = —dh ou t/dr/i at some point at the outflow 
and the vectors dh out /drji lie in the tangent plane to dVL out while b is nowhere in this plane. We have 
moreover that L 6 C^r^). Indeed, we can prove that all the derivatives of L are bounded. In order to 
do it, let us remark that the differential of X (£',£(£')) represents a vector in the tangent plane at some 
point on dfl ou t so that it is perpendicular to the outward normal n. We have thus for any i = l,...,d — 1 

= „.(g ({ ., i ( f . ))+ g (f , £ ( 5 .»|( 5 . ) ) = „.(g (f , £ ( 5 -)) +Hf ., £K0 ,|(O 

so that 

n>b(t'M?)) 

and this is bounded since X has bounded partial derivatives and n-b>a by the hypothesis. Note also 
that L is strictly positive. 

• We can now establish the decomposition V = G® ± A- Take any 0eV nC 1 (ri) and introduce 
peL 2 (Q) by 

M ^ ^ J Lia <K?,t)j(t>,t)dt 

p{x)=p((, ,£d)=p(u= tlaF) , — ; ■ L2 ) 

f Lii) j(?,t)dt 

(from now on we switch back and forth between the Cartesian coordinates x= (xi,...,x d ) and 
the new ones (£',£d) = (£1, ■ ■ • >£d))- Evidently, p is constant along each field line. Moreover, p is 
the L 2 -orthogonal projection of cf> on the space of such functions. Indeed, if ip = ip(£') 6L 2 (f2) 
is any function constant along each field line then 

P ^dx= f p^jdt= f P (ew) / m 

jd Jr in Jo 

HO r 

<K? MM?) = / Udx. 

Jn 
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Let us prove that p € V, i.e. that its derivatives are square integrable. The change of variable 
t = L(£')s yields the function 

P fij(£',L(Z')s)da 

Now we have dp/d^d = and for all 9p/9&,t=i,...,d-i denoting a = a(£') = (L J(£' ,L(£')s)ds)~ 1 , 
<p = (/>(£' ,L(£')s) and same for J we obtain 



9 P 9(1 [\jA j_ JA j. f 1 W 91 JA 

— — = —— I (pJas + a / ——Jds + a / — — ttt-s J as 



96 



« Jo 



o d£ 



96i 96 



+a / (p-^—ds + a j (f>-^r-^rsas 



o 9& 



(1.3) 



o 



9£ d % 



Using all the previous bounds on the functions L and J and skipping the details of somewhat 
tedious calculations, we arrive at 



dp 

96 



dx = 
<C 



r,„ Jo 



L ^ (dp 
U6 
HO ( 



r,„ Jo 



Jdtddd' 
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0& 



implying 



9p 
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<c 



L 2 (fl) 



lz, 2 (Q) ' 



d<j) 


2 


d<j> 




96: 


+ 

L 2 (n) 


96i 


2 ) 

L 2 (n)/ 



<cu\\ 2 m{n) 



Thus p£i/ 1 (ri), hence and q = (f> — p&A. Since the dependence of p on </> is continuous 

in the norm of H 1 ^), a density argument shows that the decomposition <f>=p-\-q with p€G 
and g G A exists for any e V. 

Let us now introduce the operator P as the L 2 -orthogonal projector on Q, that means 

P-.V^G, <j>EVi — >P<j>eg givenby (TO) . 

Then, the estimates in the preceding paragraph show that the operator P is continuous in the 
norm of i? 1 (SI) : 



||Vx(P0)||£» ( n)<C||V0|| i3(n)l V^eV 
• We have also the following Poincare-Wirtinger inequality: 

H0-^IU 2 (fi)<C||V||0|| L2(o) , V(/>e V. 



(1.4) 
(1.5) 



To prove this, it is sufficient to establish that | \q\ < C||V||<?| |.L 2 (fi) for all g £ A We observe 
that 



and 



Inll 2 

\1\\L 2 (Q) ■ 



\L 2 (fl) : 



r<„ Jo 
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(6,6i)J(C',6i)^dC- 



The requirement q £ A is equivalent to 



I q(e,Zd)J(Z',ZdHd=o f-a.a.e'er m . (1.6) 

Jo 

We have thus to prove for every £' 

provided (jl.6p . Fixing any making the change of integration variable £d = L(£')t and intro- 
ducing the functions u(t) =q(£',L(£')t)J(£' ,L(£')t) and J(t) = J(£',L(£')t), we rewrite the last 
inequality as 

Since L u(t)dt — we have by the standard Poincare inequality 

L rl 

u 2 {t)dt<C 2 P I (u'(t)) dt. (1.8) 



• Let us turn to the verification of Hypothesis B'. Take any u € V. We want to prove that one can 
decompose u—p + q with p 6 Q and q^C and the trace of u on d£li n (denoted g) is in L 2 (dQ,i n ). 
In the ^-coordinates we can write a surface element of dflm as da = S 1 with a function S 

smoothly depending on £ . We see now that for u suffuciently smooth 



\\9\\l^ n) = I 9 2 {Os(Z')d£,' 



r-l 

<c ' 



(by a one-dimensional trace inequlity) 

<c||u||£ 

By density, the trace g is thus defined for any m£V with ||s l ||z,2(gfj i „) <C|Mly- Taking p = 
p(£') = g(£') we observe by a similar calculation that ||p||l 2 (£"2) < C|Mly so tnat By 
definition q = u—p<E C. 

Appendix B. On the choice of the finite element space Cu- 

Let Q be the rectangle (0,L X ) x (Q,L y ) and the anisotropy direction be constant and aligned with 
the y-axis: b= (0, 1). Let us use the Qfc finite elements on a Cartesian grid, i.e. take some basis function 
9 Xi {x), i = 0,...,N x and 9 yj (y), j = 0,...,N V and define the complete finite element space Xh (without 
any restrictions on the boundary) as spa,n{9 Xi (x)8 y . (y) < i < N x , 0<j< N y }. The following subspace 
is then used for the approximation of the unknowns p,q,l€V 

Vh = {vhEX h /v h \au D =°}- 
We want to prove that taking for the approximation of A,/iG£ the space Ch under the form 

C h = {\ h £X h /\ h \ d n in =0}, (2.1) 
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leads to an ill posed problem (|3.1[) . 



Claim There exist AftG£ft, Aft=^0, such that a^(Xh,Ph) = for all PftGVft. In fact there are exactly 
2N y linearly independent functions having this property. 

Remark B.l. In the continuous case, the equation 



implies A = by density arguments. These density arguments are lost when discretizing the spaces V 
resp. C. 

Proof of the Claim. We can suppose that the basis functions 9ij(x,y) :=9 Xi (x)8 yj (y) are enumerated 
so that 8ij(0,y) = for all z > 1 and 0oj(0,y) ^0. Hence for all Ph = J2Pij®ij ^V/i, the coefficients satisfy 
Poj =0 since the part of the boundary {a; = 0} is in dflu- Let M — (m,ik)o<i.k<N x be the mass matrix 
in the x-direction: rriik — J 6 Xi (x)8 Xk (x)dx. This matrix is invertible, hence there is a vector a£l" I+1 
that solves Ma = e with e€R Nx+1 , e = (1,0, . . .,0)*. Take any fixed integer j, l<j<N y and define 
Aft e C h as A/,. = J2 a iQi 3 • Then, for all p h = J^Pki 9kl S V h we have 



As we can do this for all (i,j), i = 0, 1 <j<N y and in the same manner for all i = N x , 1 <j<N y , 

there are 2N y linearly independent functions with the property aii(Aft,£>ft) = for all ph £ Vft. 

We see now that the system (|3.ip with zero right hand side / = possesses non-zero solutions 
(p e h , \ e h , q\, 1^, = (0,A^, 0,0,0) where \ 3 h is any of the functions constructed in the preceding para- 
graph. It means that (|3.1[) is ill posed, i.e. the corresponding matrix is singular. 
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